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Abstract 

A  method  is  presented  for  constructing  geometric  design  data  from  noisy  3-D  sen¬ 
sor  measurements  of  physical  parts.  In  early  processing  phase,  RLTS  regression  filters 
stemming  from  robust  estimation  theory  are  used  for  separating  the  desired  part  of 
the  signal  in  contaminated  sensor  data  from  undesired  part.  Strategies  for  producing 
a  complete  3-D  data  set  from  partial  views  are  studied.  Multiple  representations  are 
used  in  model  construction  because  there  is  no  single  representation  that  would  be  most 
appropriate  in  all  situations.  In  particular ,  surface  triangulation ,  NURBS ,  and  superel¬ 
lipsoids  are  employed  in  order  to  represent  efficiently  polygonal  and  irregular  shapes , 
free  form  surfaces  and  standard  primitive  solids.  The  size  of  the  required  control  point 
mesh  for  spline  description  is  estimated  using  a  surface  characterization  process.  Sur¬ 
faces  of  arbitrary  topology  are  modeled  using  triangulation  and  trimmed  NURBS.  A 
user  given  tolerance  value  is  driving  refinement  of  the  obtained  surface  model.  The 
resulting  model  description  is  a  procedural  CAD  model  which  can  convey  structural 
information  in  addition  to  low  level  geometric  primitives.  The  model  is  translated  to 
IGES  standard  product  data  exchange  format  to  enable  data  sharing  with  other  pro¬ 
cesses  in  concurrent  engineering  environment.  Preliminary  results  on  view  registration 
using  simulated  data  are  shown.  Examples  of  model  construction  using  both  real  and 
simulated  data  are  also  given. 


1  Introduction 


In  this  paper  we  present  an  approach  for  integrating  an  intelligent  sensory  system  into 
a  part  of  design  automation  system.  Solid  modelers  could  benefit  geometric  models 
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constructed  automatically  and  rapidly  from  3-D  sensory  data.  Such  tools  are  useful  as 
a  design  aid,  especially  for  modeling  free  form  shapes  which  is  a  very  time  consuming 
design  task  by  hand,  and  requires  extensive  knowledge  about  the  modeling  tools,  such 
as  splines.  Sometimes  no  design  data  exists  for  an  old  part  and  the  redesigning  could 
be  done  by  reverse  engineering  the  part  from  sensor  measurements.  Customizing  is  also 
often  needed,  and  it  is  desirable  to  keep  the  unit  price  affordable  although  the  number  of 
parts  to  be  produced  is  small.  The  analysis  of  the  part  and  process  planning  could  also 
be  started  in  very  early  phase  of  the  design  process  using  the  initial  model  constructed 
from  the  sensor  measurements. 

The  two  main  research  problems  we  are  facing  in  the  CAD  model  construction  from 
3-D  sensor  data  are: 

1.  Data  acquisition  and  combination  of  the  partial  data  sets  into  a  complete  3-D 
data  set. 

2.  Data  interpretation  by  fitting  models. 

The  first  problem  requires  estimation  of  relative  rotation  and  translation  between  data 
set  obtained  from  different  vantage  points  and  combination  of  all  data  into  one  common 
coordinate  frame.  The  goad  of  data  interpretation  is  to  produce  a  geometric  model  of  a 
part  to  be  imported  into  a  solid  modeling  system.  Similarly  to  Computer  Aided  Geo¬ 
metric  Design  (CAGD),  there  is  no  single  method  or  representation  in  Computer  Vision 
that  would  be  appropriate  in  all  situations.  Therefore,  we  employ  multiple  represen¬ 
tations  in  model  construction.  The  produced  geometric  model  should  be  compatible 
with  common  representations  in  modeling  systems  in  order  to  analyze  and  simulate 
the  model  and  share  it  with  other  automation  subsystems.  The  designer  should  also 
be  able  to  modify  the  model  because  the  design  typically  evolves. 

Our  approach  constructs  procedural  CAD  models ,  which  are  procedures  that  gen¬ 
erate  the  part  geometry.  Procedural  models  are  able  to  represent  low  level  geometry 
of  the  part  as  well  as  its  overall  structure.  Structural  information  is  vital  for  analysis, 
simulation  and  process  planning  and  it  must  be  detected  by  these  processes  if  not  pro¬ 
vided  by  the  geometric  model.  Moreover,  procedural  models  are  useful  in  representing 
intersections  of  surfaces,  for  example,  in  the  case  of  trimmed  parametric  surfaces.  The 
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intersection  is  described  in  the  procedure  and  it  can  be  approximated  in  the  level  of 
required  accuracy  when  it  is  actually  needed.  The  designer  is  also  able  to  modify  the 
procedure,  if  necessary. 

The  capability  to  communicate  between  different  subsystems  during  the  design 
process  is  a  prerequisite  for  concurrent  engineering.  The  data  sharing  is  provided  by 
standard  product  data  formats,  such  as  IGES  [5,  12].  The  proposed  system  is  depicted 
as  a  part  of  concurrent  engineering  environment  in  Figure  1. 


CAPP 

\  / 


Figure  1:  The  proposed  system  as  a  part  of  a  concurrent  engineering  en¬ 
vironment.  The  CAX  processes  are  Design,  Engineering,  Process  Planning, 
Manufacturing  and  Inspection. 

The  organization  of  this  paper  is  as  follows.  In  section  2  we  address  data  acquisition, 
view  registration  and  integration  problems.  In  section  3  we  describe  briefly  shape 
representations  used  in  model  construction.  In  section  4  we  show  some  examples  using 
real  and  simulated  range  data.  Finally  in  section  5  we  summarize  and  discuss  some 
areas  requiring  future  research. 


2  Data  acquisition 


We  chose  to  use  optical  non-contact  sensors  for  measuring  3-D  shape  of  the  objects. 
Coordinate  Measuring  Machines  (CMM)  were  not  considered  because  of  their  low  speed 
in  acquiring  data  from  free  form  shapes  which  require  dense  measurements. 

Physical  measurements  are  prone  to  errors.  In  the  case  of  range  data,  the  actual 
noise  distribution  differs  from  the  nominal  one  which  is  often  Gaussian.  Furthermore, 
there  mav  occur  outliers  due  to  the  orientation  of  material  of  the  surface  or  because 


of  other  statistical  populations  present  in  the  processing  window.  The  raw  data  is 
filtered  using  RLTS  [13,  14]  robust  regression  filters  in  order  to  recover  the  structure  of 
the  underlying  signal  and  reject  outliers  which  may  cause  incorrect  estimates  in  model 
building  processes. 

In  general,  optical  non-contact  active  range  data  acquisition  techniques  provide 
incomplete  3-D  information  because  the  signal  does  not  reach  all  the  surface  points  if 
the  data  is  obtained  from  one  viewpoint  at  the  time.  A  complete  3-D  data  set  has  to 
be  merged  from  a  collection  of  images  from  different  viewpoints.  The  rotation  groups 
of  regular  polyhedra,  as  noted  in  [4],  provide  a  convenient  set  for  uniformly  sampling 
the  observation  sphere.  Therefore,  the  scanning  procedure  should  use  such  a  set  of 
evenly  distributed  viewpoints  as  a  default,  if  no  symmetry  is  obvious. 

In  order  to  combine  multiple  range  views  into  one  complete  3-D  data  set  the  regis¬ 
tration ,  i.e.,  the  relative  rotation  and  translation  between  the  views,  must  be  estimated, 
and  the  integration  of  the  views  into  nonredundant  data  set  in  a  common  coordinate 
frame  performed.  Recent  overviews  of  the  research  on  view  registration  and  integration 
methods  are  given  in  [30,  4]. 

2.1  Registration  and  integration 

2.1.1  Background 

The  registration  estimates  the  relative  transformations  between  different  views  and 
transforms  all  the  data  into  a  common  coordinate  frame.  Typically,  methods  assume 
either  that  the  transformation  is  known,  or  corresponding  features  are  detected  reliably 
from  each  view  and  subsequently  the  transformation  can  be  solved  accurately.  In  the 
latter  case,  the  features  are  a  set  of  a  priori  known  reference  points  from  the  environ¬ 
ment  that  are  visible  in  different  views,  or  features  extracted  from  object  surfaces.  This 
approach  is  adequate  in  simple  situations  where  the  object  consists  of  relatively  few 
geometric  primitives  that  can  reliably  detected  from  different  viewpoints.  In  the  case 
of  sculptured  free-form  shapes,  however,  it  is  difficult  to  establish  correspondencies. 

We  chose  to  adapt  the  Iterative  Closest  Point  (ICP)  algorithm  proposed  by  Besl  and 
McKay  [4]  for  registering  a  single  view  with  a  known,  computer-generated  database. 
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It  was  chosen  because  no  feature-to-feature  correspondencies  are  required,  it  is  com- 
putationnaily  efficient  and  independent  from  data  representation,  as  long  as  a  method 
for  computing  point/primitive  distances  exists.  Its  main  shortcoming  is  obviously  that 
only  a  locally  optimal  displacement  is  found. 

The  method  matches  a  collection  of  points  from  one  set  of  raw  data  with  a  set  of 
primitives  from  a  model.  Each  point  is  basically  associated  to  its  closest  primitive, 
the  type  of  primitive  defining  the  exact  distance  measuring  function.  Given  a  model 
M=  {mk}  containing  M  primitives,  let  V=  {pi}  be  the  set  of  N  points.  We’ll  see  how 
to  choose  those  points  in  the  next  section.  X=  {ii}  is  then  defined  as  the  “projection” 
of  V:  Xi  is  the  closest  point  of  the  closest  primitive  mkt  to  pi.  Once  a  match  V  — ►  X  is 
established,  an  optimal  displacement,  in  the  form  of  a  registration  state  vector  q,  which 
consists  of  a  rotation  quaternion  qn  and  a  translation  vector  qj,  can  be  computed.  The 
mean  square  error  function  [4]: 

^  1=1 

is  minimized.  Let’s  call  qo  the  state  vector  solution,  qo  is  applied  on  the  set  V  and 
the  closest  primitive  search  is  iterated  on  the  displaced  set,  leading  to  a  new  V  — +  X 
matching.  The  consecutive  displacements  q  obtained  are  cumulated  into  a  final  solution 
until  the  error  measure  e(q)  is  less  than  a  certain  threshold. 

2.1.2  Constrained  solution 

The  view  registration  problem  in  model  construction  context  is  more  complicated  be¬ 
cause  there  is  no  a  priori  known  part  model.  Hence,  we  have  to  incrementally  build  the 
internal  representation  used  as  the  part  model,  by  adding  each  view  after  registration 
instead.  On  the  other  hand,  we  obtain  an  approximate  transformation  between  the 
views  because  the  movement  of  sensors  or  objects  are  controllable.  Such  approximate 
transformation  can  serve  as  a  good  first  estimate  for  the  ICP  procedure  in  order  to 
find  a  global  minima  instead  of  a  local  one.  Furthermore,  it  provides  a  way  to  exclude 
points  from  the  matching  process  if  they  are  not  visible  in  the  representation  registered 
so  far. 

The  main  problem  is  that,  using  a  given  set  Vk  extracted  from  range  image  k ,  some 
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points  pi  may  not  appear  in  the  model  M  compiled  so  far  (from  images  0  — »  k  —  1). 
In  that  case,  some  matches  in  the  correspondence  V  — ►  X  will  always  be  erroneous, 
leading  to  a  wrong  estimation  of  q  as  pointed  out  in  [4], 

Fortunately,  if  a  reasonable  estimate  of  q  is  known  prior  to  the  beginning  of  the 
registration  procedure,  overlapping  points  will  produce  relatively  small  errors,  whereas 
mismatched  points  often  lead  to  large  residual  errors,  and  can  then  be  detected  by 
analyzing  the  distribution  of  the  errors: 

II Xi  -  R{qR)pt  -  9r||2, 

where  q  is  the  state  vector  applied  to  the  data  in  the  previous  iteration.  To  further 
separate  outliers  from  valid  points,  we  check  for  the  consistency  of  the  normals  after 
displacement: 

•  Rq-R(nft)  >  cos(a ) 

where  n. $  is  the  normal  estimated  at  point  u,  and  a  a  threshold  angle.  The  normals 
can  be  estimated  by  using  local  window  operators  or  robust  estimation  techniques  if 
the  image  is  noisy  [13].  Non-valid  points  display  large  consistency  errors  and  can  be 
discarded  at  each  iteration:  At  iteration  A;,  the  displacement  q  is  recomputed  until  no 
point  violates  the  normal  consistency  constraint,  a  is  kept  large  to  maintain  flexibility 
in  the  matching  process  and  insure  convergence.  Also  it  is  clear  that  such  a  method 
will  fail  if  gross  outliers  (non  overlapping  data)  are  predominant  in  the  X  set  and  the 
first  estimation  of  q  prior  to  refinement  is  wrong. 

Another  adaptation  concerns  the  choice  of  points  and  primitives.  Many  types  of 
primitives  are  possible,  depending  on  the  geometry  of  the  scanned  surface.  Prelimi¬ 
nary  tests  have  been  performed  successfully  using  large  planar  patches.  We  are  now 
investigating  the  case  of  free  form  objects  using  Delaunay  triangulated  [11]  surface 
representation  addressed  later  in  shape  representation  section.  To  run  the  ICP-based 
algorithm  on  triangulated  data,  we  select  the  set  V  as  the  centers  of  gravity  of  the 
triangles  in  one  image  (a  low  resolution  is  chosen,  so  that  typically  less  than  a  hundred 
points  are  selected).  Those  points  are  registered  against  a  triangulation  of  the  exist¬ 
ing  model,  constructed  with  the  previous  registered  images,  the  first  image  being  the 
starting  point.  Eventually  a  complete  3D  data  set  will  be  incrementally  constructed. 
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A  more  thorough  description  of  the  constrained  registration  is  given  in  [31]. 

2.1.3  Integration 

The  integration  part  combines  registered  and  partly  overlapping  data  sets  into  a  com¬ 
plete  nonredundant  3-D  data  set.  The  methods  typically  perform  a  (weighted)  av¬ 
eraging  of  the  samples  that  are  in  the  overlapping  parts  of  the  views.  The  level  of 
integration  and  representation  of  the  data  have  to  be  chosen  as  well.  The  approaches 
are  either  data  driven,  such  as  [29]  where  the  surfaces  are  represented  using  low  level 
primitives,  e.g.,  triangular  mesh,  or  model  driven,  where  more  elaborate  model  sur¬ 
faces  are  employed  and  the  surface  parameters  are  adjusted  to  fit  to  the  data.  An 
example  of  latter  approach  is  given  in  [6],  where  initial  triangulated  ellipsoid  or  icosa¬ 
hedron  is  refined  to  approximate  the  bounding  surface.  Other  possible  approaches  are 
described,  for  example  in  [19,  27].  In  our  case,  a  data  driven  integration  part  is  under 
development. 

3  Issues  on  shape  representation 

The  type  of  shape  representation  have  to  be  chosen  based  on  not  only  the  family  of 
shapes  we  are  describing  but  also  on  the  task  where  the  representation  is  used.  In  object 
recognition  from  single  arbitrary  viewpoint,  properties  such  as  viewpoint  invariance  and 
uniqueness  are  important.  In  geometric  modeling,  however,  the  representation  have  to 
be  unambiguous  but  not  necessary  unique  and  the  geometry  is  typically  described  in 
object  centered  coordinate  frame. 

Constructive  Solid  Geometry  (CSG)  and  Boundary  Representation  (B-rep)  are 
widely  used  in  CAGD  [23].  CSG  models  standard  primitive  solids  effectively  but  mod¬ 
eling  of  sculptured  free  form  surfaces  is  difficult.  B-rep  defines  a  solid  by  its  bounding 
surfaces.  Polygonal  representation  is  commonly  used  for  modeling  flat  surfaces.  Non 
Uniform  Rational  B-splines  (NURBS)  are  widely  used  for  modeling  free  form  shapes, 
and  they  are  also  able  to  represent  conics  and  quadric  surfaces  exactly.  Sweep  meth¬ 
ods  are  often  used  to  design  solids  that  have  rotational  or  translational  symmetry  [23]. 
The  ’’Design  by  Features”  paradigm  [26]  uses  manufacturing  features  to  create  the 
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part  geometry  which  is  very  appealing  approach  because  the  design  is  often  done  for 
manufacturing.  It  seems,  however,  that  there  is  no  single  representation  or  method 
that  would  be  best  for  every  design  task.  Therefore,  several  CAD-systems  are  hybrid 
systems,  i.e.,  use  multiple  representations  to  provide  efficient  tools  for  different  design 
tasks  [23]. 

The  representations  used  for  CAD-based  computer  vision  can  be  generally  classi¬ 
fied  into  volumetric,  surface  and  sweep  representations  [3].  Volumetric  methods  are 
actually  surface  based  methods  because  surface  evaluation  is  required  to  recover  vol¬ 
umetric  description.  Such  methods  can  represent  only  closed  surfaces,  hence  they  are 
not  suitable  for  describing  partial  information.  Surface  based  methods  represent  the 
part  geometry  typically  by  a  set  of  bounding  low  order  surface  patches.  Generalized 
cylinders  (GC)  are  a  typical  sweep  representation.  It  seems,  as  in  the  design,  that  there 
is  no  single  representation  that  would  recover  an  appropriate  description  from  sensor 
data  in  all  situations. 

In  order  to  model  different  shapes  we  are  employing  multiple  representations.  An 
optimal  triangulation  is  generated  for  modeling  polygonal  and  irregular  shapes  and 
surfaces  of  arbitrary  topology.  Triangulation  includes  very  little  structural  information 
about  the  part  but  it  can  be  used  as  a  worst  case  representation,  if  no  other  method  is 
suitable.  NURBS  are  used  for  modeling  free  form  surfaces  because  of  their  continuity 
and  local  control  properties.  Furthermore,  they  are  included  in  IGES  product  data 
exchange  standard  which  facilitates  data  sharing  and  concurrent  engineering.  Superel¬ 
lipsoid  models  are  used  to  detect  overall  part  structure  which  allows  us  to  use  more 
efficient  model  primitives  that  are  helpful  in  part  analysis  and  process  planning. 

3.1  Triangulation 

In  order  to  describe  polygonal  and  irregular  surfaces  and  surfaces  of  arbitrary  topology 
we  construct  a  collection  of  triangles  (triangulation)  describing  the  surface.  We  require 
that  the  distance  between  any  point  of  the  object  surface  and  its  projection  onto 
the  triangulation  along  the  surface  normal  is  less  than  a  user  defined  tolerance  value. 
Triangulation  methods  are  widely  used  in  approximation  theory,  finite  element  analysis, 
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and  CAD  [27]. 

A  k/n-triangulation  is  a  collection  of  A;-simplices  in  n- dimensional  space.  A  triangle 
is  a  2-simplex  and  a  tetrahedron  is  a  3-simplex,  for  example.  In  general,  triangulations 
where  triangles  are  nearly  equilateral  are  desirable.  In  particular,  we  are  generating  a 
Delaunay  triangulation  which  satisfies  the  property  that  the  interior  of  the  minsphere 
of  a  A;-simplex  contains  no  vertex  of  any  fc-simplex.  The  minsphere  is  the  smallest 
(n  -  l)-dimensional  sphere  which  passes  through  k+  1  vertices  of  a  ^-simplex.  In  two- 
dimensional  case,  the  minsphere  property  means  that  the  minimum  angle  is  maximal 
over  all  triangulations. 

The  method  applied  here  follows  the  Generalized  Delaunay  Triangulation  procedure 
presented  in  [11]  with  k  =  2  and  n  =  3.  The  method  starts  with  one  fc-simplex  and 
new  triangulation  points  are  inserted  one  at  the  time.  Different  insertion  operations 
are  executed  based  on  the  location  of  insertion  point  relative  to  the  triangulation.  The 
process  continues  until  all  points  have  been  inserted  or  or  the  remaining  points  cannot 
be  inserted  without  violating  the  empty  minsphere  property. 

The  triangulation  is  refined  based  on  Euclidean  distances  between  points  and  tri¬ 
angles  by  adding  and  removing  points  to  meet  the  user  defined  tolerance  value  tol. 
For  point  removal  as  well  as  point  addition,  all  is  needed  is  to  compute  the  distance 
between  a  point  and  a  triangle  in  3D.  A  point  is  candidate  for  addition  if  the  distance 
d  to  the  nearest  triangle  T  is  >  tol.  The  list  of  points  to  be  added  is  obtained  by 
sampling  the  input  data  at  various  resolution,  starting  at  a  coarse  level,  and  interating 
until  no  point  can  be  added.  Between  each  resolution,  the  triangulation  vertex  list  is 
scanned  for  possible  removals:  For  each  point  p  in  the  triangulation,  let  T  be  the  set 
of  2-simplices  sharing  p .  p  is  tentatively  removed,  and  T  is  retriangulated  into  T* .  p 
is  definitively  removed  only  if  the  distance  of  p  to  T7  is  lesser  than  tol. 

Given  a  triangle  T  =  (pi,P2,P3),  the  distance  between  a  point  p  and  primitive  T  is 
given  by: 

d(p,  T)  =  min  ||aipl  +  a2p2  +  a3p3  -  p\\ 

ai+a2+a3=l 

with  the  weighting  coefficients  d\  6  [0,1],  a<i  £  [0,1],  a3  £  [0,1].  In  practice,  one 
projects  p  on  the  plane  embedding  T,  and  the  projection  pi  should  be  inside  T.  When 
the  triangulation  is  generated,  one  can  easily  compute  the  equation  of  the  edges  of  T. 
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Let’s  then  call  d,(v)  the  algebraic  distance  between  any  point  v  and  edge  i  ( i  is  the  edge 
which  doesn’t  contain  p]).  To  perform  the  interiority  check,  one  has  only  to  compute 
d,(pl)  (*  €  {1,2,3}),  then  tests: 

0  <  di(p1 )  <  d,(pi)  for  i  €  {1,2,3} 

provided  that  d,(pi)  >  0  (otherwise,  reverse  the  sign  of  the  inequalities).  The  d,(p*) 
are  already  computed  and  stored  during  the  triangulation  process. 


3.2  Spline  approximation 


B-spline  surfaces  have  several  desirable  properties  for  geometric  modeling  and  high 
quality  surface  approximation,  such  as  local  control  and  continuity.  B-spline  surfaces 
lie  within  the  convex  hull  formed  from  the  control  point  mesh.  It  is  important  to  have 
sufficient  number  of  control  points  to  be  able  to  describe  all  the  degrees  of  freedom 
of  the  underlying  surface.  If  there  are  too  few  control  points  the  fit  is  not  likely  to 
converge  and  on  the  other  hand,  the  fewest  number  of  control  points  tends  to  yield 
the  fairest  surfaces  [25].  An  initial  estimate  for  the  B-spline  control  point  mesh  size 
is  computed  by  using  the  maximum  number  of  geometrically  homogeneous  patches  in 
each  parameter  direction.  The  surface  patches  are  detected  by  a  local  characterization 
process.  If  the  order  of  B-splines  is  three  (degree=2),  we  need  at  least  three  control 
points  in  each  parameter  direction  to  be  able  to  describe  each  second  order  surface 
patch. 

A  non-uniform  rational  B-spline  surface  (NURBS)  is  a  more  general  case  of  nonra- 
tional  B-spline  surface,  and  is  defined  as  a  bivariate  polynomial  function  of  parameters 


«  and  r  as  follows: 


S(u,v)  = 


£?=i  £7=1  h^N^{u)Mhl{v) 


'  ’  '  £"=1  £7=1  ’  v  ' 

where  Nhk  and  M}j  are  the  basis  functions,  h{j  are  the  weights,  and  the  B,,j’s  are  the 
control  points,  n  and  m  identify  the  number  of  control  point  vertices  in  each  direction. 
The  basis  functions  Nt^  of  order  k  are  defined  recursively  as  follows: 


NiM 


1  ,if  Xj  <  u  <  x)+i 
0  ,  otherwise 
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■Nitk(u)  = 


(3) 


(u- Xi)Ni<k-i(u)  (xj+k  -  u)Ni+hk-i(u) 

%i+k— 1  ”  3'i  %i+k  —  ^i+l 

where  xt’s  are  ordered  set  of  knots  from  knot  vector.  A  convention  0/0  =  0  is  used 
for  the  basis  function  computation.  Basis  functions  Mjj  of  order  l  for  parameter  v  are 
computed  similarly.  We  chose  to  employ  maximum  of  4th  order  (k  =  l  =  4,  cubic) 
B-spLines  to  make  the  approximation  less  sensitive  to  small  local  variations. 

Chord  length  parameterization  is  employed.  The  parameter  values  are  normalized 
to  [0, 1]  range.  An  open  end  condition  is  used  to  force  the  spline  to  begin  exactly  from 
the  first  control  point  and  end  at  the  last  control  point. 

The  locations  of  the  control  points  of  the  approximating  B-spline  surface  are  com¬ 
puted  by  minimizing  errors  in  least  squares  sense.  Now  we  have  to  solve  ’s  from 
equation  (1),  and  5(u,  v)Js  are  the  measured  data  points.  All  the  weights  are  originally 
set  to  1.0  because  5(u,  v)’s  are  physical  measurements.  The  weights  of  the  control 
points  can  be  adjusted  later  in  surface  refinement  [20].  Using  matrix  representation 
the  solution  is: 

[B]  =  [[Cf  [CT'lCf  [S].  (4) 

where  elements  of  C  are  Cij  =  S  is  the  matrix  of  data  points,  and  B  is  the 

obtained  control  point  mesh. 

The  accuracy  of  the  approximation  should  meet  a  user  given  tolerance  value.  The 
error  of  the  approximation  is  defined  as  the  Euclidean  distance  between  the  measured 
and  the  approximated  surface  with  same  (u,v)  parameter  values.  There  are  different 
approaches  for  controlling  the  shape  the  B-spline  surface,  see  [8,  17,  20,  25].  We 
start  with  a  good  estimate  of  the  appropriate  control  point  mesh  size  and  add  knots 
if  the  distance  exceeds  the  given  tolerance  value.  Adding  a  certain  number  of  knots 
has  a  consequence  of  adding  the  same  number  of  control  points.  Curve  or  surface 
discontinuities  can  be  introduced  by  inserting  a  knot  with  multiplicity  equal  to  the 
order  of  the  B-spline. 

Tensor  product  B-spline  surfaces  require  a  rectangular  parametric  grid  which  maps 
the  coordinates  in  (u,v)  parameter  space  to  physical  (x,y,z)  coordinates.  Resampling 
may  be  necessary  to  get  a  rectangular  arrangement  of  the  data.  In  the  case  of  scattered 
samples,  for  instance,  the  points  are  typically  organized  into  triangular  faces  and  then 
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the  contours,  e.g.,  isoparametric  lines,  are  interpolated  [19]. 

It  is  not  possible  to  describe  certain  surfaces  of  arbitrary  topology  with  single 
nondegenerate  B-spline.  Designers  can,  however,  can  introduce  degeneracies  into  the 
mesh  by  reducing  the  control  points  forming  an  edge  of  the  mesh  into  a  single  point  or 
use  non-tensor  product  patches  [16].  The  surfaces  can  then  be  represented,  for  example, 
by  using  Gregory  patches  or  Rational  Boundary  Gregory  (RBG)  patches  [7]  which  are 
joined  together.  The  continuity  properties  at  the  joints  are  relaxed  into  geometric 
tangent  plane  (G1)  continuity.  In  order  to  transfer  such  data  to  CAD  systems  using 
IGES  [12]  product  data  exchange  format  the  representation  have  to  be  converted  to 
NURBS. 

A  simple  engineering  solution  for  situations  where  the  rectangular  arrangement  of 
tensor  product  surfaces  is  not  appropriate  is  to  employ  trimmed  surfaces.  A  trimmed  B- 
spline  surface  is  essentially  a  regular  B-spline  surface  where  certain  parts  of  the  surface 
are  marked  ’’invalid”  [10].  In  our  approach,  the  boundaries  of  the  surface  are  used  to 
compute  trimming  curves  which  divide  it  into  invalid  and  valid  parts.  The  intersection 
curves  of  parametric  surfaces  are  computed  using  subdivision  based  techniques  [8].  The 
fit  procedure  is  run  using  a  bounding  box  for  the  object,  and  the  parts  of  the  surface 
which  are  not  on  the  object  surface  are  declared  invalid.  Trimmed  surfaces  are  included 
in  IGES  standard  and  advanced  solid  modelers  [12,  1]. 


3.3  Superellipsoid  model  recovery 


Superellipsoids  are  a  subclass  of  superquadrics  [2]  that  can  represent  shapes  ranging 
from  ellipsoids  to  cuboids  and  cylindroids.  A  superellipsoid  surface  in  nonparametric 
implicit  form  is  defined  as  follows  [28]: 


f(x,y,z) 


(5) 


where  aj,  a2,  and  a3  define  the  size  in  x-,  y-  and  z-axis  direction.  £\  and  e2  are  the 
shape  (squareness)  parameters  in  the  latitude  and  in  the  longitude  plane,  respectively. 
Additional  parameters  are  employed  to  describe  global  deformations.  The  Levenberg- 
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Marquartd  method  [22]  is  used  for  minimizing  the  following  expression  [28]: 

_N_ 

rnm^{y/a^ai(F(xi,yi,Zi)  -  l)]2,  (6) 

»'=i 

where  the  function  F(x,y,z )  =  f(x,y,z)£l  determines  the  locus  of  a  point  relative 
to  superellipsoid  surface  and  N  is  the  number  of  samples.  The  range  of  the  shape 
parameters  is  constrained  to  0.1  <  £i,£2  <  2.0.  The  goodness  of  fit  is  defined  as 
GOF  =  v/a1a2a3((E^:i  2/>  z)  ~ 

The  superellipsoid  model  recovery  is  used  to  detect  typical  primitive  solids  shapes 
and  overall  part  properties  such  as  symmetry.  Rotationally  symmetric  shapes  are  con¬ 
structed  using  surface  of  revolution  design  primitive  whereas  translationally  symmetric 
shapes  can  be  generated  by  extrusion.  The  primitive  solids  [18]  and  the  approximating 
superellipsoid  shape  parameters  are  depicted  in  Table  1.  The  nature  of  the  obtained 

Table  1:  Primitive  solids,  the  corresponding  CAD  model  parameters,  and  the 
approximating  superellipsoid  shape  parameters.  means  that  in  general  the 
primitive  can  not  be  recovered  using  the  superellipsoid  model  we  employ. 


superellipsoid  parameters  is  qualitative  and  they  can  be  used  as  a  hypothesis  to  invoke 
the  appropriate  model  building  procedure.  If  the  shape  parameters  indicate  that  the 
part  is  a  natural  quadric,  a  more  accurate  description  is  obtained  by  fitting  quadric 
model  to  the  data  [21].  In  the  case  of  a  surface  of  revolution,  the  actual  model  building 
process  fits  conic  sections  along  the  assumed  axis  of  symmetry  to  recover  the  rotation 
axis  and  the  profile  NURBS  curve  accurately. 


Primitive  solid 

sphere, ellipsoid 

parallelepiped 

cylinder 

cone 

torus 

wedge 


Model  parameters 

radius,  major  and  minor  axis 

length,  width,  height 

height,  radius,  normal  plane 

base  circle,  vertex  point 

major  and  minor  radius,  initial  plane 

length,  width,  height 

length,  width,  height 

vertex  points 

regularized  intersection  of  surfaces 


Shape  parameters /remarks 

5j  =  1.0,  t2  =  1.0 

e1  <  1.0,52  <  1.0  or  5j  =  2.0,52  <  1.0 

5!  -  1.0,52  <  1.0 _ 

cylinder  +  tapering 
super  tori 

only  a  subset,  parellelepiped  +  tapering 
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4  Experimental  results 


In  the  experimental  part  we  show  preliminary  results  on  view  registration  using  simu¬ 
lated  data.  Moreover,  model  construction  results  are  shown  using  simulated  and  real 
sensor  data  emphasizing  the  need  for  multiple  representations. 

The  images  are  filtered  using  robust  RLTS  filtering  to  recover  the  structure  of  the 
original  signal  and  reject  outliers.  Very  deviating  observations  would  cause  serious 
errors  in  model  construction  which  employs  least  squares  error  norm  in  fitting.  Figure 
2  shows  filtering  results  using  a  5  point  processing  window  for  a  simulated  sample  profile 
where  Gaussian  noise  with  p  =  0  and  a  =  5.0  and  random  bit  error  with  probability 
P  —  0.015  are  added  to  the  noise-free  signal.  A  12  bit  quantization  is  used. 


Figure  2:  Filtering  results  for  noisy  signal:  (a)  The  noisy  signal,  and  (b)  the 
RLTS  filtered  signal,  respectively. 

Registration  estimates  the  relative  rotation  and  translation  between  the  views.  In 
our  preliminary  experiments,  the  ICP  algorithm  is  used  to  minimize  the  distances  of 
the  points  to  planar  surface  patches.  A  region-growing  algorithm  uses  surface  normal 
consistency  to  generate  a  planar  approximation  of  the  part.  Small  patches  that  occur 
typically  in  the  vicinity  of  Co  and  C\  discontinuities  are  discarded  as  unstable  and  irrel¬ 
evant.  The  centers  of  gravity  of  the  remaining  planar  areas  are  taken  as  reference  points 
Pi,  to  be  matched  against  the  planar  primitives  extracted  in  a  contiguous  viewpoint, 
using  the  ICP  based  registration  scheme.  A  simulation  using  a  computer-generated 
polyhedral  object  in  eight  different  poses  is  shown  in  Figure  3. 

The  final  registration  is  shown  with  cross-section  along  the  x,y  and  z  axis  of  the 
object  (Figures  4  and  5).  The  algorithm  was  here  able  to  recover  the  right  displacement 
in  all  cases  with  good  accuracy:  1/4  pixel  accuracy  in  pixel  to  pixel  registration  between 
views.  Apparently,  less  accurate  results  are  obtained  using  real  data.  Furthermore,  one 
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has  to  take  account  the  actual  amount  of  overlap  between  views.  The  more  complicated 
the  data,  the  more  occlusions,  and  the  higher  the  number  of  viewpoints  is  required  to 
obtain  a  complete  data  set. 


15 


Figure  4:  The  registration  result  is  illustrated  by  taking  a  cross-section  of  the 
object  by  keeping  X  coordinate  (=  1.3)  constant.  At  the  top  is  the  complete 
data  set,  the  other  images  represent  the  contributions  of  each  view.  The  white 
area  is  due  to  the  cutting  plane  grazing  of  one  of  the  object  surfaces. 
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Figure  5:  A  collection  of  isocontours  along  the  X  axis  (1.05  <  X  <  1.5).  For 
each  isocontour,  data  were  collected  from  all  the  views  using  the  registrations 
obtained  via  the  ICP  algorithm.  The  first  view  was  used  as  a  reference  frame. 
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The  model  construction  examples  are  given  using  the  data  sets  depicted  in  Figure 
6.  The  Cylindrical  Pin  data  is  synthetic  data  where  Gaussian  distributed  noise  with 
zero  mean  and  a  =  4.0  as  well  as  random  bit  error  noise  with  probability  P  =  0.001  are 
added.  The  image  resolution  is  256-by-256  pixels  and  16  bit  quantization  is  used  for 
the  depth  values.  The  Face  Mask  image  and  the  Hand  image  from  NRCC  [24]  range 
image  library  are  chosen  to  demonstrate  the  capability  to  model  free-form  surfaces  and 
surfaces  of  arbitrary  topology. 


Figure  6:  Test  images:  a)  The  Cylindrical  Pin  is  a  synthetic  range  image 
produced  from  a  CAD  model,  b)  the  Face  Mask  and  c)  the  Hand  are  real 
range  images  from  NRCC  range  image  library. 

The  surface  triangulation  is  produced  in  order  to  model  surfaces  of  arbitrary  topol¬ 
ogy,  polygonal  shapes,  and  irregular  shapes  that  do  not  consists  of  smooth  surfaces 
typically  found  in  manufactured  objects.  Furthermore,  such  representation  is  useful 
in  view  registration  and  integration.  Triangulation  does  not  convey  much  structural 
information  about  the  part  geometry.  The  starting  point  of  the  triangulation  was  cho¬ 
sen  on  the  center  of  gravity.  The  triangulations  were  refined  by  adding  or  deleting 
points  until  a  tolerance  value  was  achieved.  In  the  case  of  the  Face  mask,  724  triangles 
were  needed  to  meet  the  0.5  mm  tolerance  value  for  the  description.  The  Hand  data 
required  1260  triangles  with  tolerance  value  0.4  mm.  Triangulations  of  the  Face  Mask 
and  Hand  are  depicted  in  Figure  7. 

The  superellipsoid  model  recovery  is  used  to  reveal  global  shape  properties.  The 
obtained  shape  parameters  are  used  as  a  hypothesis  to  invoke  the  appropriate  model 
building  procedure.  Table  2  shows  the  recovered  shape  parameters  and  the  quality 
of  the  fit  measures  for  the  test  images.  The  shape  parameters  reveal  the  rotational 
symmetry  of  the  Cylindrical  Pin.  The  quality  of  the  fit  is  also  high,  hence  surface  of 
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Figure  7:  Triangulation  examples:  a)  The  Face  Mask,  and  b)  the  Hand. 


revolution  modeling  primitive  is  selected.  The  quality  of  the  fit  is  low  for  the  Face 
Mask  and  the  Hand  data  and  they  are  modeled  as  a  collection  of  bounding  surfaces. 
In  addition,  a  tapering  deformation  is  obtained  for  the  Hand  data.  The  recovered 


Table  2:  Superquadric  model  recovery  results 


Parameters/ 

Test  image 

Shape 

Goodness  of  fit 

GOF 

Cylindrical  Pin 

Cl  =  0.10,52  =  108 

0.06 

Face  Mask 

c\  =  0.15,52  =  0.87 

0.15 

Hand 

5 1  =  0.10,  52  =  0.10 

0.22 

superellipsoid  models  for  the  test  pieces  are  depicted  in  Figure  8. 


Figure  8:  The  obtained  superellipsoid  models  of  the  test  pieces:  a)  the  Cylin¬ 
drical  Pin,  b)  the  Face  mask,  and  c)  the  Hand 

Sculptured  surfaces  are  approximated  using  NURBS.  Surface  characterization  is 
employed  to  estimate  the  number  of  control  points  needed  in  the  control  mesh  to  be 


able  to  describe  all  the  degrees  of  freedom  of  the  underlying  surface  using  B-splines.  For 
each  second  order  surface  patch  we  need  3  control  points  in  each  parameter  direction. 
The  maximum  number  of  control  points  in  each  parameter  direction  was  selected  to  be 
the  size  of  the  surface  mesh  in  that  direction.  A  patch  is  considered  significant,  hence 
included  in  the  mesh  size  estimation,  if  the  number  of  pixels  exceeds  a  given  threshold 
value.  Surface  characterization  results  for  test  images  are  depicted  in  Figure  9. 


Figure  9:  Surface  characterization  results:  a)  the  Cylindrical  Pin,  b)  the  Face 
Mask,  and  c)  the  Hand. 

We  chose  to  use  chord  length  parameterization  in  spline  approximation.  In  general, 
the  accuracy  is  good  if  the  surface  is  smooth.  Larger  errors  are  caused  by  rapid  changes 
in  the  surface  shape.  The  surface  description  is  refined  to  meet  a  user  defined  tolerance 
value  by  inserting  knots,  and  consequently  control  points.  A  knot  is  added  to  a  point 
where  local  error  maxima  exceeding  the  tolerance  value  occurs.  Figure  10  depicts  a 
profile  from  the  Face  Mask  and  the  corresponding  error  distances  to  the  approximating 
spline  obtained  using  chord  length  parameterization.  The  result  of  refinement  by  knot 
insertion  using  tolerance  value  LO  mm  are  also  shown. 

Surface  discontinuities  cause  errors  because  of  the  continuity  property  of  B-splines 
are  violated.  Weaker  continuity  properties  can  be  introduced  by  inserting  multiple 
knots  at  same  point,  for  instance,  four  in  the  case  of  subdividing  a  cubic  B-spline 
where  a  discontinuity  occurs.  An  example  of  subdivision  is  depicted  in  Figure  11. 

Sometimes  rectangular  arrangement  of  data  is  completely  inappropriate.  For  such 
situations  we  chose  to  use  trimmed  surfaces.  The  fit  is  run  on  the  bounding  box 
of  the  part  and  the  boundaries  are  used  for  computing  trimming  curves.  A  fairly 
dense  control  point  mesh  is  needed  to  isolate  the  errors  introduced  on  the  boundaries. 
Two  different  parameterization  methods  were  experimented  in  the  context  of  trimmed 
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Figure  10:  a)  The  original  profile  from  Face  Mask  and  the  obtained  B-spline 
(dashed  line),  and  b)  the  corresponding  error  distances,  c)  shows  the  B-spline 
refinement  result  using  1.0mm  tolerance  value,  and  d)  the  error  distances, 
respectively. 

surfaces:  Chord  length  parameterization  consumes  most  of  the  parameter  space  by 
the  boundaries,  whereas  uniform  parameterization  provides  less  accurate  result  by 
the  boundaries  but  it  gives  better  approximation  for  the  valid  parts  of  the  surface. 
Therefore,  use  uniform  parameterization  and  insert  more  knots  where  large  errors  occur 
in  order  to  make  the  approximation  more  accurate.  An  example  of  surface  trimming 
operation  is  performed  on  the  Hand  image.  The  boundaries  are  detected  using  Deriche 
edge  detector  [9].  The  complete  B-spline  surface  and  the  resulting  trimmed  surface  are 
depicted  in  Figure  12. 


Figure  11:  B-spline  subdivision  allows  representing  discontinuities:  The  sim¬ 
ulated  data  from  the  Cylindrical  Pin  and  the  obtained  B-spline  (dashed  line) 
a)  before  and  b)  after  the  subdivision. 
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Figure  13:  The  CAD  model  for  a)  the  Cylindrical  Pin,  and  b)  and  c)  for  the 
Face  Mask,  respectively. 
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Figure  14:  Model  data  for  the  Cylindrical  Pin:  a)  a  part  of  the  Alpha.l  model 


procedure, 


and  b)  a  part  of  the  IGES 


description. 
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5  Summary 


In  this  paper,  we  presented  a  computer-aided  engineering  tool  for  building  geometric 
models  of  parts  from  3-D  sensor  data.  Data  have  to  be  acquired  from  several  viewpoints 
and  integrated  into  a  complete  3-D  data  set  in  common  coordinate  frame.  We  showed 
preliminary  results  on  registration  using  a  simple  polyhedron  type  object.  Multiple 
representations  are  used  in  model  construction  in  order  to  model  efficiently  different 
shapes  and  consequently,  employ  appropriate  CAD  modeling  primitives.  Moreover, 
there  is  no  single  representation  that  would  be  always  appropriate.  In  particular,  we 
choose  to  employ  superellipsoids,  NURBS  and  Delaunay  triangulation  in  order  to  cover 
standard  geometric  shapes  as  well  as  free  form  and  irregular,  complicated  shapes.  The 
result  is  a  procedural  CAD  model  which  is  able  to  convey  structural  information  about 
the  part.  A  procedure  which  is  generating  the  part  geometry  is  relatively  easy  to 
modify  by  the  user  which  is  necessary  because  the  design  typically  evolves.  The  model 
building  is  addressed  as  a  part  of  concurrent  engineering  environment.  Hence,  the 
model  have  to  be  translated  to  standard  product  data  exchange  format  to  enable  data 
sharing  with  other  processes. 

The  ongoing  and  future  research  is  directed  toward  refining  and  extending  the  data 
acquisition  part  in  order  to  combine  data  from  different  vantage  points  even  if  the  object 
has  sculptured  free  form  shape.  In  particular,  surface  triangulations  are  employed  to 
register  the  views  incrementally.  Methods  for  modeling  surfaces  of  arbitrary  topology 
are  also  developed  further.  Moreover,  interfacing  to  engineering  analysis  tools,  e.g., 
dynamic  simulation,  is  under  work.  The  motivation  is  to  be  able  to  view  a  part, 
accurately  model  and  represent  it,  and  predict  the  effect  of  different  design  parameters 
on  the  performance  of  the  mechanical  system  in  order  to  optimize  or  redesign  the  parts 
without  actually  building  prototypes. 
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